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O ■ Abstract 
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The dissipative quantum dynamics of an anharmonic oscillator coupled to a bath is studied with 
\ the purpose of elucidating the differences between the relaxation to a spin bath and to a harmonic 

bath. Converged results are obtained for the spin bath by the Surrogate Hamiltonian approach. 
■ This method is based on constructing a system-bath Hamiltonian, with a finite but large number 

> : . 

t^J- , of spin bath modes, that mimics exactly a bath with an infinite number of modes for a finite time 

interval. Convergence with respect to the number of simultaneous excitations of bath modes can 

<n : 

be checked. The results are compared to calculations that include a finite number of harmonic 
modes carried out by using the multi-configuration time-dependent Hartree method of Nest and 
Q_i! Meyer, [J. Chem. Phys. 119, 24 (2003)]. In the weak coupling regime, at zero temperature and 

for small excitations of the primary system, both methods converge to the Markovian limit. When 
initially the primary system is significantly excited, the spin bath can saturate restricting the energy 
| acceptance. An interaction term between bath modes that spreads the excitation eliminates the 

•i-H . 

saturation. The loss of phase between two cat states has been analyzed and the results for the spin 

S ■ 

and harmonic baths are almost identical. For stronger couplings, the dynamics induced by the two 
types of baths deviate. The accumulation and degree of entanglement between the bath modes 
have been characterized. Only in the spin bath the dynamics generate entanglement between the 
bath modes. 



*Electronic address: davg@fh.huji.ac.il 
^Electronic address: christiane.koch@lac.u-psud.fr 
^Electronic address: ronnie@fh.huji.ac.il 



1 



I. INTRODUCTION 



Modeling quantum many-body systems is a challenging problem. The main obstacle is 
the exponential growth in complexity with the number of degrees of freedom. Significant 
simplifications are achieved by partitioning the total system into a primary part and a 
bath describing the environment [if. The idea is to model the primary system explicitly 
and the bath implicitly, thus minimizing the complexity of the bath to its influence on the 
primary system. A bath composed of a set of noninteracting harmonic oscillators is the one 
most widely used. The idea originates from a normal mode analysis combined with a weak 

n 

system-bath coupling assumption |2| . If the bath is only weakly perturbed by the system, 
it can be considered linear, and therefore described as a collection of harmonic oscillators. 
Such a bath is natural for systems interacting with the radiation field [3|. The harmonic 
bath model has also been applied to less favorable scenarios such as energy relaxation and 
dephasing of molecules in the liquid phase or on solids. In these cases a strong coupling or 
interactions with a low-temperature environment may cause large system-bath correlations, 
and will therefore result in a failure of the Markovian approximation. To overcome such 
difficulties in the dynamics of molecules that are in intimate interaction with an environment, 

n 

an alternative approach termed the Surrogate Hamiltonian |4j has been developed. The 
Surrogate Hamiltonian method employs a bath composed of two-level systems that acts as 
a spin bath Q E , Q, B, Q] 

The concept of the system-bath separation underlines the quantum description of many 
body dynamics. The origins of the spin and harmonic baths are different. The harmonic bath 
is closely related to a normal mode decomposition. Once this is done the spectral density 
function is able to completely determine the relaxation dynamics. From a computational 
point of view the determination of the spectral density is a major task. The most popular 
working procedure is to extract it from classical mechanics ^(|. The drawback is that this 
procedure assumes harmonic modes and a linear system bath coupling term. The spin bath 
has its origin in a tight binding model of condensed phase. This can also become a simulation 

Hcedure if the parameters of the tight binding model can be estimated from first principles 

The purpose of the present study is to compare the performance of the two baths in 
a simple system composed of a primary anharmonic oscillator coupled to a multi-mode 
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bath. In the limit of weak system-bath coupling, it has been shown that the two baths are 
equivalent. For finite temperature the equivalence requires a rescaling of the spectral density 
'unction which determines the coupling of the primary system to the different bath modes 
3 0, Q]- The limiting coupling strength where the dynamics induced by the two baths 
differ has not yet been characterized. For stronger coupling strength, the ergodic behavior 
of the two baths should be different. The bath modes of the linearly driven harmonic bath 
are uncorrelated. In the spin bath the coupling to the primary system induces quantum 
entanglement between the different modes. It is valuable to know how this fundamental 
difference influences the dynamics of the primary system. 

Our comparative study is based on a numerical model of a system coupled to a bath, with 
a large but finite number of modes. For a finite interval of time determined by the inverse 
of the energy level spacing, the finite bath mimics exactly a bath with an infinite number of 
modes. For this interval the primary system cannot resolve the full density of states of the 
bath. By renormalizing the system-bath interaction term to the density of states, the finite 
bath faithfully represents the infinite bath up to this time limit. 

The dynamics of the primary oscillator coupled to the harmonic bath has been re- 
cently calculated based on the multi-configuration time dependent Hartree approximation 



(MCTDH) 12,ll3j. The authors were able to show that for a Morse oscillator coupled to a 



bath, converged results could be obtained for a bath consisting of 60 modes to a time scale 
of 3 ps. The present study utilized the same system and system-bath coupling parameters, 
but employed a spin bath in the context of the Surrogate Hamiltonian. The comparison 
allows an evaluation of the similarities and differences between the two descriptions. Once 
the differences are identified, it becomes possible to modify the Surrogate Hamiltonian bath 
to extend the realm of similarity. 

3oth cases is not Markovian and differs from the Redfield 
ll7 , Q| . In the weak coupling limit the numerical study 



The system-bath construction in 



16 



14 , 13 or semigroup treatments 

n 

of Nest and Meyer |12| was able to identify a coupling parameter where the Markovian 
semigroup limit was reached. One can reason that the Surrogate Hamiltonian bath should 
behave similarly in this range of coupling parameters. 

The present paper is organized as follows: Section |H] outlines the theory of the two 
models: the Surrogate Hamiltonian approach which employs a bath of two-level systems 
(TLS) and the MCTDH method using an harmonic bath. Section UTTl describes the system 
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used for calculations. Section IIVI compares the results for the two different environments. 
The standard process investigated in studies of quantum dissipative dynamics is energy 
relaxation (cf. Section TlV Ajl . An interaction between bath modes is introduced in Section 
IIVBI The difference between correlated and uncorrelated initial states is the subject of 
Section IIVCI In addition, the decoherence in the TLS bath of the Surrogate Hamiltonian 
is compared to that in a bath of harmonic oscillators (cf. Section Ft V D|) . A characterization 
of different kinds of entanglement in the Surrogate Hamiltonian approach is presented in 
Section HV El Finally, Section IVl summarizes and concludes. 

Appendix [X] compares the equations of motion between the two different types of bath 
and Appendix El introduces two different measures of entanglement of a two-spin system. 
It should be noted that atomic units are used throughout the paper (h = m e = a = 1) . 

II. THEORY 

The system under study describes a primary system immersed in a bath. The state of 
the combined system-bath is described by the wave function \I/(R, /3 1; . . . ,/3 2 N ) where R 
represents the nuclear configuration of the dynamical system, and {/3j} are the bath degrees 
of freedom. The Hamiltonian of such a combined system is: 



H = H s <g> I B + Is ® H B + H SB 



(1) 



The primary system Hamiltonian takes the form: 



H s = T + y s (R) , 



(2) 



where T = P /2M is the kinetic energy and Vs is an external potential, which is a function 
of the system coordinate(s) R. Hb denotes the bath Hamiltonian consisting of an infinite 
sum of single mode Hamiltonians h\,: 




j 



For the harmonic bath the single mode Hamiltonians take the form: 




(4) 
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where Pj,Qj are the normal mode momentum and coordinate respectively, and a.j = 
\ I 3 !? 1 q, H , 1 p.- is the corresponding annihilation operator. For the spin bath: 

V 1 J y/2m,jWj ■> 

hj = ujj&l&j , (5) 

where <xL &j are the standard spin creation and annihilation operators of mode j. 

The system-bath interaction Hsb can be decomposed into a sum of products of system 
and bath operators without loss of generality. Specifically a system-bath coupling inducing 
vibrational relaxation is considered: 

H S b = -/(R)®^V j , (6) 

j 

where Vj = Xfy = Aj(aj + a,-) for the harmonic bath and Vj = Aj(<rj + <Xj) for the spin 
bath. /(R) is a function of the system coordinate operator. The influence of the bath on 
the primary system is characterized by the spectral density function J(u). To include the 
density of states, the definition of the spectral density function is chosen as 0, [3] : 

J (u) = ^2\\j\ 2 p(u)5(u -uJj) , (7) 
j 

that is the system-bath coupling is weighted by density of states. Thus the constants Aj are 
determined as: 

Xj = yjj{uj)/p{uj) , (8) 

where piyjj) = (ujj + i — ujj)^ 1 is the density of the states of the bath. 

Observables associated with operators of the primary system are determined from the 
reduced system density operator: 

p s (i?, R') = tiB {\^ )( *&\}, where tr^{ } is a partial trace over the bath degrees of freedom. 
The system density operator is constructed from the total system-bath wave function and 
only this function is propagated. 

Since within a finite interval of time, the system cannot resolve the full density of bath 
states, it is sufficient to replace the bath modes by a finite set. The sampling density in 
energy of this set is determined by the inverse of the time interval. The finite bath of iV 
spins is constructed with a system-bath coupling term, which in the limit iV — > oo converges 
to the given spectral density of the full bath. The Surrogate Hamiltonian, as well as the 
MCTDH method, consist of a finite number of bath modes, and they are therefore limited 



to representing the dynamics of the investigated system for a finite time (shorter than the 
Poincare period at which recurrences appear [2^). These recurrences are caused by the 
finite size of the bath so that after some time the energy flow into the bath is reflected at 
its boundaries. 

The Surrogate Hamiltonian contains all possible correlations between the primary system 
and the environment. The combined system-bath state is described by a 2 N dimensional 
spinor with iV being the number of bath modes. The spinor is bit ordered, i.e., the jth bit 
set in the spinor index corresponds to the jth TLS mode, which is excited if the counting 
of bits starts at j = 0. The dimension 2 N results from the total number of possibilities to 
combine two states N times. Thus the total wave function can be written as 



(N\ /JV\ 



|^(R,{/3 j })) = c o |0 o (R)) + ^c J |^(R))+ ^c jfc |^ fc (R)) + ... , (9) 

j=0 j,k=0 

where |0j(R)) = (0, . . . , 0j(R), . . . , 0) T is a singly-excited spinor, |0 jfc (R)) = 
(0, . . . , 0j(R), • • • , 0fc(R), • • • , 0) T is a doubly-excited spinor and so on. The jth compo- 
nent corresponds to the jth TLS being excited. However, considering all 2 N possibilities of 
combining the bath modes might not be necessary in a weak coupling limit. In this case, for 
short time dynamics, it is possible to restrict the number of simultaneous bath excitations 



2l|. As an extreme example, only single excitations might be considered. If one restricts 
the number of simultaneous excitations, the dimension of the spinor becomes the sum of 
binomial coefficients X]/£o (^zTO w ^h iV exc the number of simultaneous excitations. The 
construction is similar to the configuration-interaction (CI) approach in electronic structure 
theories. The restriction of simultaneously allowed excitations leads to significant numerical 
savings and its validity can be checked by increasing iV exc . 

In the MCTDH method (2^] the wavefunction which describes the dynamics of a 
system with M degrees of freedom, is expanded as a linear combination of time-dependent 
Hartree products: 

ni nM M 

l*(Qi, • • • ,Q*,t)) = E • • • E A n,-,nAt) II \<Pj?iQ*>t)) > ( 10 ) 

31=1 j M =l K=l 



where \<p[: ) is the single-particle function (spf) for the k degree of freedom and the Aj u 



■,3 m 



denote the MCTDH expansion coefficients. The total number of coefficients Aj lr .. jM and 
basis function combinations scales exponentially with the number of degrees of freedom M. 



Considering a system coupled to a multi-mode bath, the use of the multicon fiqurational wave 
function ensures the correct treatment of the system-bath correlations I23I. I^. The method 
also enables grouping of several modes together, which reduces both the number of single- 
particle degrees of freedom and the correlation effects between different modes. Although 
the exact treatment is contained in the limit of an infinite number of configurations, in the 
weak coupling limit, the time-dependent basis employed in the MCTDH method should 
be relatively small. Worth et al. [_23| have pointed out that even for weak coupling, one 
spf per bath mode (the Hartree limit) is not sufficient to fully describe the system-bath 
interaction. However, the number of spf 's for the bath degrees of freedom can be increased 
until convergence is achieved, which makes this approximation controllable. 



III. THE MODEL 



The primary system is constructed from an anharmonic (Morse) oscillator of mass M: 

~ 2 

H S = — + D (e"** - 2e-^) . (11) 

The coupling term is non-linear in the Morse oscillator coordinate R, but reduces to a linear 
one for a small R: 

1 _ p-oiR 

/(R) = • (12) 

a 

The spectral density function was chosen to be the same as in the harmonic bath case. For 
an Ohmic bath the damping rate 7 is frequency-independent and the spectral density in the 
continuum limit is given by 

J{uj) = M^u (13) 

for all frequencies u> up to the cutoff frequency u c . A finite bath with equally spaced sampling 
of the energy range was used. 

The parameters used are the same as in Ref. a well depth D of 0.018 a.u., a = 2 

a.u., and a mass of M = 10 5 a.u. The initial state was chosen to be a Gaussian displaced 
by Rq = 2R from the origin with a width of o = R (R « 0.09129 a.u. is the characteristic 
length scale of the Morse oscillator). For such a displacement the coupling term ()12|) is 
almost linear. The initial system-bath state has a direct product form where the bath is at 
zero temperature. Such a state has no initial correlations between the system and the bath. 
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There are a few characteristic time scales of the system. The period of the Morse oscillator 



is r osc = 27r/f2 m 127 fs, where f2 = ot\j2D/M refers to the harmonic frequency of the 
potential. The bath has two time scales. Tb a th is associated with the highest frequency 
uj c = 2.5Q and corresponds to a time scale of 52 fs. The time scale corresponding to the 
frequency spacing Au defines the Poincare period (r rec ). It should be larger than any other 
time scale of interest. With u c fixed this time becomes: 

2tt 2txN 

Tree ~ ~T = • 14 

Au oj c 

Thus, with an increasing number of bath modes, the convergence progresses in time. In 
our simulations the number of TLS is chosen to be iV = 20 ... 60 (for different coupling 
strengths), which ensures that r rec is greater than the overall simulation time. 

The calculations were performed in three different interaction regimes identified by con- 
sidering the involved time scales: (i) weak coupling referring to 7 _1 = 1630 fs ^> r osc ,rb at h; 
(ii) the intermediate situation characterized by 7" 1 = 163 fs ~ r osc > Tbath; (hi) the strong 
coupling regime defined by 7 _1 = 54 fs ~ r bath < r osc . 

In the simulations discussed below, the average position of the oscillator and the energy 
relaxation were calculated for all three coupling strengths. For comparison, the effective 
subsystem energy was defined as in jl^ |: 

E s = (H s ) = (H s ) +0.5(H SB >. (15) 

It includes half of the system-bath interaction term. 

The dynamics of the system combined with the bath is generated by solving the time- 
dependent Schrodinger equation: 

vI/(R,{/3 j },t) = e- jfe vI/(R,{/3 j },0). (16) 

Each spinor component ipj (R) is represented on a spatial grid. The kinetic energy operator 



is applied in Fourier space employing FFT 
to compute the evolution operatoi 
already been given in Ref. Q, El- 



and the Chebychev method 26] is used 
to oo-npute the evotton orator. N„, detai, s of applying the hath operators have 
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FIG. 1: The energy relaxation (lower panel) and damped oscillations of the average position (upper panel) 
of the Morse oscillator in the weak coupling limit (7 _1 = 1630 fs). The bath is assumed to be Ohmic 
with cutoff frequency oj c = 2.9 • 10~ 3 a.u. and consists of N — 60 TLS. The initial state was chosen to 
be a Gaussian displaced by Rq — 2R with a width of a = R, where R w 0.09129. Thick solid lines refer 



Q) 



to MCTDH calculations with a bath of harmonic oscillators (adopted from Ref. [12J). Dashed lines refer 
to Surrogate Hamiltonian calculations with only single excitations. Thin lines refer to two simultaneous 
excitations allowed. 



IV. RESULTS AND DISCUSSION 

A. Energy relaxation and small amplitude motion. 

First a restricted Surrogate Hamiltonian is applied, which limits the possible system-bath 
correlations. The most extreme restriction includes only single excitations. The results for 
the weak coupling case (7 _1 = 1630 fs) are shown in Fig. For a short period of time the 
energy relaxes with the same rate in the two types of bath. However, after t s m 500 fs the 
rate decreases and eventually the system energy becomes constant. It should be pointed out, 



9 




-0.0160 



si -0.0165 

H 

e -0.0170 



-0.0175 







300 600 
time (fs) 



900 





1 1 






MCTDH 






— N=40(l) 






N=40(2) 










\v ***• 


















1,1, 



900 



FIG. 2: The energy relaxation (lower panel) and damped oscillations of the average position (upper panel) 
of the Morse oscillator in the intermediate coupling regime (7 _1 = 163 fs). The bath parameters and the 
initial state are the same as in the weak coupling calculations. The number of bath modes is N = 40. 



.12]) 



Thick solid lines refer to MCTDH calculations with a bath of harmonic oscillators (adopted from Ref. 12] ). 
Dashed lines refer to Surrogate Hamiltonian calculations with single excitations only. Thin lines refer to 
calculations with two allowed simultaneous excitations. 

that the saturation time is not the recurrence (Poincare) time (t s < r rec ). This is confirmed 
by the fact that for time t > t s the overall energy transfer from the bath back to the system is 
not complete. Calculating the population of the bath modes shows that at t > t s most of the 
system energy is transferred to very few (or even one) bath modes, which are in resonance 
with the system's frequency. Modes which are near to the resonance mode or modes become 
saturated and start to transfer the excitation back to the system. A dynamic "steady state" 
between the system and the bath is formed, where most of the modes transfer energy back, 
while one (or very few) continue to absorb energy from the system. 

When the number of simultaneous excitations is increased to two, the effect of saturation 
appears at a later stage (t s > 2000 fs). The results become similar to those of Ref. [12] and 
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time (fs) time (fs) 

FIG. 3: (Left panel) The relative difference in the effective subsystem energy ((i?s) + 0.5(.Hsb)) between 
the bath with only single excitations allowed {E\) and the bath with two simultaneous excitations (-E2) as a 
function of time. The difference is calculated for a few coupling strengths. The simulation time is t = 900 fs 
and the number of bath modes is N = 10 . . . 40. (Right panel) The population P s j m of 1 and 2 simultaneous 
bath excitations is compared to the bath with two simultaneous excitations. The solid lines refer to the 
weak coupling limit (7 _1 = 1630 fs) and the dashed lines refer to medium coupling (7 _1 = 163 fs). 

the values of the average position (see Fig. Q (upper panel)) are nearly indistinguishable. We 
conclude that for the weak coupling case, the bath that has two simultaneous excitations is 
completely sufficient to reproduce the dynamics generated by all simultaneous excitations 
for times up to 2 ps. 

The relaxation dynamics for medium coupling are shown in Fig.|2j A saturation effect was 
obtained for the bath restricted to single excitations. However, two simultaneous excitations 
were sufficient to overcome this saturation and converge the whole dynamics of the problem. 
A slight difference in the energy relaxation rate of the two baths is identified. The TLS 
bath causes stronger relaxation, but the results are still in good agreement with those of 
Ref. |l2j. Since the initial state is a function of R and the system-bath coupling depends on 
R as well, the initial excitation influences the effective strength of the coupling. If the initial 
displacement, i.e. the initial excitation of the primary system, is decreased, the saturation is 
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FIG. 4: The energy relaxation (lower panel) and damped oscillations of the average position (upper panel) 
of the Morse oscillator in the strong coupling strength (7 _1 = 54 fs). The bath parameters and the initial 
state are the same as in the weak coupling calculations. The number of bath modes is N = 20. Solid lines 



refer to Nest and Meyer's MCTDH calculations (adopted from Ref. 12]), where the thick lines are the 
full-dimensional wave packet result and the thin lines refer to the Morse-Lindblad model. Dashed lines refer 
to Surrogate Hamiltonian calculations with single excitations. Dashed-dotted lines refer to a bath with five 
excitations allowed simultaneously. 

postponed. We can then deduce that the relaxation rate converges to the value of Ref. [12J • 
Combining the results of Figs. ^ and 121 leads to the conclusion that the differences between 
the two types of bath in the weak and intermediate coupling regimes are caused by the 
saturation of a few "central" modes in the spin bath. This saturation is postponed if the 
bath includes more correlations. For very weak coupling, these higher order system-bath 
correlations become insignificant. 

The problem of including all system-bath correlations is therefore crucial in the medium 
and strong coupling regime. Fig. El shows the difference in the system energy ((Hg) + 
0.5 (H SB » for two bath with only single excitations and a bath in which two simul- 
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taneous excitations are allowed. The calculations were made for different coupling strengths. 
As the coupling strength is reduced, the difference decreases. Thus in a very weak coupling 
limit, the TLS bath with only single excitations (no system-bath correlations) becomes suffi- 
cient to describe the dynamics for relatively long times. In this limit the TLS bath coincides 
completely with the harmonic bath. 

The issue of including system-bath correlations has also been addressed in the MCTDH 
calculations. In Ref. the same system has been studied with the G-MCTDH method (the 
MCTDH with Gaussian expansion functions). Differences between the single- configur at ional 
(the Hartree limit) and the multi-configurational descriptions (with an increasing number 
of single particle functions) have been obtained for the energy relaxation process. In these 
calculations at least four single particle functions per resonant bath modes and two spf for 
secondary modes were required to achieve convergence in the relatively weak coupling limit 
( 7 -i = 500 fs). 

In the strong coupling regime (Fig. 0} there is considerable deviation between the two 
models. In the Surrogate Hamiltonian model the energy relaxes faster and the oscillator 
is damped after a single period. The relaxation rate obtained for the TLS bath is closer 



to the one obtained by the Morse-Lindblad model (adapted from Ref. 12]). As expected 
convergence requires many simultaneous excitations. For example, a bath consisting of 
N = 20 modes required at least five simultaneous excitations for converging the energy 
relaxation dynamics. 

B. The interaction between the bath modes. 

The saturation of the TLS bath modes can be eliminated by allowing energy exchange 
between bath modes. This is done by adding to the bath Hamiltonian Hb of Eq.Q the 
term: 

Hint = ^ K V^\^3 + ' ( 17 ) 

where the parameter «y(= is the interaction strength between two bath modes. The 
interaction can be restricted to the nearest neighbors in energy by the condition Kij = for 
|i — j\ > 1. The detailed algorithm of applying Eq.([17|) in the bit representation has been 
described in Ref. |27[ in the context of pure dephasing. 
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The term cr\&j + o")o"t describes a two quasi-particle interaction within the bath. A 
qualitative picture is based on an almost elastic exchange of energy between the two nearest 
neighbor bath modes which are almost degenerate. The process is described by a creation 
of an excitation in one mode at the expense of another and vice versa. 

The new bath Hamiltonian including the interactions can be diagonalized leading to: 

H B = y^a^g-j , (18) 

i 

where cu, are the eigenvalues of D^(Hb + Hi nt )D. In the new basis of {cri} the system-bath 
interaction term in Eq. (jUJ) is also modified. However for sufficiently small k the eigenvalues 
of the bath change only slightly, but the saturation effect is postponed to a much later time. 

Fig. El shows the influence of the interaction between the bath modes on energy relaxation 
in the weak coupling limit = 1630 fs). The dynamics are calculated for a relatively 
long period of 3 ps. For such a long time the saturation effect is observed even for a bath 
with two simultaneous excitations. Since the saturation time t s is determined mostly by the 
saturation of the few modes close to resonance with the subsystem, increasing the number 
of modes cannot prolong t s . 

Adding an interaction between the bath modes leads to slower decay and delayed satu- 
ration (Cf Fig. 03 upper panel). This can be understood from the following considerations: 
the interaction term, Eq. (jl7|) . describes the transport of excitation from one bath mode to 
its nearest neighbor. Consequently, k determines how quickly the excitation is transported 
away from a TLS mode close to resonance with the primary system. On the other hand, the 
interaction energy, i.e. the expectation value of (Hsb), depends on the population of the 
primary system and of the bath modes close to it. If the population is removed from those 
bath modes and "diffuses" all over the bath, the interaction energy decreases and the decay 
becomes slower. This explains the upper panel of Fig. [5] which shows the energy relaxation 
for different values of k. 

An optimal value of n can minimize the differences between the spin and the harmonic 
bath. To demonstrate the effect, the calculations were carried out for iV = 60 bath modes 
and k — 1.5 • 10~ 4 (Cf. the lower panel of Fig. EJ). For this value of k the spectrum of the 
bath was only slightly altered (less than 1%). The energy relaxation in this case is almost 



indistinguishable from the results obtained in Ref. 
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FIG. 5: The energy relaxation with interaction between the bath modes is shown in the weak coupling limit 
(7 _1 = 1630 fs). The bath is assumed to be Ohmic with cutoff frequency ui c — 2.9 • 10~ 3 a.u. (Upper panel) 
The influence of the parameter k is shown for a bath of N = 40 modes. (Lower panel) The energy relaxation 
is shown with the optimal parameter of k = 1.5 • 10 -4 (thin line). The thick solid line refers to MCTDH 
calculations with a bath of harmonic oscillators (|l2j). The dashed line refers to Surrogate Hamiltonian 
calculations with two simultaneous excitations allowed without interaction between the bath modes. The 
bath consists of N = 60 modes. 

C. Correlated versus uncorrelated states. 



ly uncorrelated system-bath state is not consis- 



The widely used assumption of an initia 
tent with most experimental situations [sgl^T 30|. The influence of initial correlations has 
been addressed in the context of the weak coupling approximation, where it appears as an 
additional inhomogeneous term [sflj ]. 

A fully correlated initial state is easily obtained in the Surrogate Hamiltonian method. 
Once the system-bath Hamiltonian H is set, the correlated ground state can be determined 
by propagating an initial guess wave function in imaginary time using H j^lj]. 
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time (fs) system-bath coupling r| 



FIG. 6: Effect of initial correlations. The energy relaxation and damped oscillations of the average position 
are shown for the initially uncorrelated (solid lines) and correlated (dashed lines) states. The dynamics 
for weak (j^ 1 = 1630 fs) and strong couplings (7 _1 = 54 fs) are compared. A bath consisting of N = 10 
modes with 2 and 5 simultaneous excitations (for the weak and strong coupling, respectively) is sufficient to 
obtain the converged results. (Left) The effective subsystem energy (upper panel) (E s = (Hs) + 0.5(-£Tsb)) 
and the expectation value for the Morse coordinate (lower panel) are shown. (Right) The bare subsystem 
energy (upper panel) E$ = (Hg) in the strong coupling regime (7 -1 = 54 fs) is shown for the short time 
dynamics. The energy stored in the system-bath coupling (lower panel) is calculated as a function of the 
coupling constant rj = M"f. 

The influence of initial correlations is shown in Fig. |H1 The uncorrelated state is identical 
to that of the previous calculations: the primary system is defined as a shifted Gaussian wave 
packet, while the bath is not excited. For the correlated initial state, the ground state of the 
total system was calculated first. Then this ground state was displaced by the shift operator 
in momentum space D = e _li?ok with Rq = 2R. The dynamics of the correlated state are 
compared to that of the uncorrelated state for weak and strong couplings (7 _1 = 1630 fs 
and 7 _1 = 54 fs). The dashed and solid lines in Fig. El correspond to the initial state being 
correlated and uncorrelated, respectively. 

The short-time dynamics differ for the correlated and uncorrelated cases, since the corre- 
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3- 



lations need to be built up in the uncorrelated case [32| . In the latter case an initial slippage 
in the system energy can be observed before the reduced dynamics appear to be Markovian 
(right upper panel). This effect is insignificant for weak coupling. Even for strong cou- 
pling, the differences between the correlated and uncorrelated cases were found to be very 
small. Apparently, the displacement is a stronger "perturbation" than that caused by the 
correlations, i.e. the displacement establishes a new initial state [33]. 

D. Decoherence. 

Decoherence has become a popular term used to describe loss of phase in coherent su- 
perpositions of quantum states due to interaction with a bath. It is therefore natural to 
compare the decoherence properties of the spin bath to those of the harmonic bath. The 
first difficulty is that there are different approaches to the definition of decoherence. Alicki 
34I ] identifies pure decoherence (dephasing) with the decay of the off-diagonal elements of 
the density operator, which is not accompanied by dissipation. He then argues that de- 
phasing cannot be caused by a harmonic oscillator bath with a coupling, which is linear in 
coordinates or momenta. 

Energy relaxation is also accompanied by loss of phase. For comparison, we will consider 
decoherence as a process caused by energy relaxation, which is characterized by a time 
T\ = 7 _1 . The decoherence effect will be illustrated in terms of the dissipative dynamics 
of cat states, defined as a superposition of two coherent states. The interaction with the 



environment leads to decay of the coherences of such a superposition on an extreme 



y short 



time scale, usually much shorter than the corresponding relaxation time scale |35j. This 
process has been modeled using G-MCTDH by |l3| . and as a result can be used for a 
comparative study. 

The Wigner function of the cat state consists of two Gaussians centered at (±-Ro,Po) and 
an interference term, which is centered at the origin. The off-diagonal part of the density 
matrix in the coherent-state basis, which contains information about quantum interferences 
between the two components of the cat state, decays with the rate 7 co h- in the Markovian 
limit the decay rate is proportional to the square of the distance between the coherent states. 
For zero temperature it is given by jsfi. I^: 

7co h = • (19) 
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FIG. 7: The decoherence effect in terms of decay of the coherence norm and off-diagonal elements in the 
energy representation. The coherence norm n co h is shown as a function of time for two different couplings, 
7 _1 = 1630 fs and 7 _1 = 500 fs. In the golden rule limit, off-diagonal elements of the reduced system density 
matrix p(R, R ) decay exponentially with the decoherence rates 7^ = 130 fs and 7^ = 40 fs (for the two 
given couplings, respectively). This decay is shown by the full lines. The dashed lines refer to the calculated 
decay of the decoherence norm. The dotted lines show the decay of trs{p^ oh } in the energy representation. 

ujq and 5 are parameters of the primary system (ujq represents the frequency of the harmonic 
oscillator, and 5 is the separation distance between the coherent states). 

The decoherence rate for a primary system coupled to a TLS bath is calculated and 
compared to the calculations for a bath of harmonic oscillators with the G-MCTDH method 
[l3l |. The calculations are performed for a cat state in a harmonic oscillator potential with 
ujq = 10~ 3 a.u. and M = 10 5 a.u. The bath has the same parameters as for the previous 
calculations with damping rates of 7 _1 = 1630 fs and 7 _1 = 500 fs. 

A quantitative measure of decoherence of the primary system is the coherence norm used 
by Strunz et al. 

n coh (t)=tT S {p c s°\t)p c s oh \t)] , (20) 

where Pg oh refers to off-diagonal elements of the subsystem reduced density matrix in the 
basis of coherent states. 

Since decoherence is a basis-dependent phenomenon, one can ask if it can also be measured 
in the basis of the eigenstates of the system Hamiltonian Hs- The question arises, whether 
these eigenstates form a pointer basis - the basis with respect to which off-diagonal 
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elements in the reduced density operator disappear due to decoherence. To perform this test 
the system density operator p s (t) = trg{p}, which has been calculated in the coordinate 
basis is transformed to the basis of the Hs eigenstates. Decomposing such a state to a 
dynamical and a static part leads to 0|: 

PsW = PcohW + C 2 Ps q (21) 

where Pg q is the equilibrium stationary system density operator and C 2 is an overlap func- 
tional given by 

C 2 = tT S {p s (t)-p^}/tT S {pf} . (22) 

Pcoh(^) m Eq. (|21j) . has no diagonal elements in the energy representation and is therefore 
traceless. Thus the decoherence effect is measured by the decay of tr s {pl oh }- 
Fig. [3 shows the decay of the coherence norm n co h and tr^jp^} for two different coupling 
strengths (both are weak). The thick lines refer to a simple exponential decay predicted 
by Eq. (jl9j) for a harmonic bath and confirmed by Ref. ^|. The dashed lines refer to the 
calculated decay of n coh , which is in good agreement with the prediction. The decay of 
trs{p^ oh } (the dotted lines) has almost the same rate at a relatively short time. However, 
the off-diagonal elements of p s in the energy representation do not decay strictly to zero. 
Therefore, in this case, the system energy eigenstates cannot be considered as a pointer basis 



E. Entanglement. 

Entanglement between two quantum states is a manifestation of additional quantum 
correlation. For example entanglement between the system and the bath means p ^ p s ®p B . 



In a dissipative environment it is expected t 
system are lost leading to decoherence [ul. 



rat initial entanglement between parts of the 



42j. In addition a bath can also provide an 



indirect interaction between totally decoupled parts of the primary system and entangle 
them Q,Q- 

The difference between the harmonic and the spin baths should be manifested in another 
type of entanglement - quantum correlations between different bath modes. A system in- 
teracting with the spin bath, can induce entanglement between two spin modes, which are 
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FIG. 8: Measurement of entanglement between the bath modes as a function of time. (Upper panel) 
The smallest eigenvalue of the partial transposition p Tj of the reduced density matrix for pair of the bath 
modes is calculated according to Appendix [5] The negative eigenvalues are averaged over all possible 
combinations of the bath modes. (Lower panel) The relative number of entangled pairs of the bath modes 
as a function of time. The calculations are performed for three different system-bath coupling strengths 
(7 _1 = 1630, 500, 163 fs). The bath consisting of N — 40 modes with two simultaneous excitations allowed 
is used in all calculations. 

not directly interacting with each other. In the harmonic bath on the other hand, a system 
linearly coupled to different modes is not able to entangle those modes (see Appendix IX]). 



Peres 



and Horodecki et al. 



46| have provided a criterion, based on partial transposi- 



tion, to determine whether a given mixed state of two subsystems is entangled (cf Appendix 
lB|) . Since the criterion is defined only for two coupled TLS, the study of entanglement is 
limited to two bath modes (i) and (j). The density operator of any two bath modes py is 
obtained as a partial trace of trk^ij{p B } over the rest of the bath modes, where the density 
operator of the bath is p B = tr^jp}. The procedure checks whether the partial transposition 
of py with respect to one of the modes has negative eigenvalues. The smallest eigenvalue 
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FIG. 9: The entanglement of formation E(p) is calculated for a couple of the bath modes The average 

over all possible pairs is shown as a function of time. The calculations are performed for three different 
system-bath coupling strengths (j^ 1 = 1630,500, 163 fs). The bath consisting of N = 40 modes with two 
simultaneous excitations allowed is used in all calculations. 

Ao of the partial transpose matrix p Tj constitutes the criteria. Then the eigenvalues with 
Ao < are averaged over all entangled pairs of bath modes. 

In Fig. |H1 the averaged parameter Ao is shown as a function of time for three different 
coupling strengths. The entanglement calculations were based on converged results obtained 
for a bath of N = 40 modes. This was sufficient to a time scale of 900 fs, for all three system 
bath coupling strength considered. Since at t — the bath is not excited, there are no 
entangled bath modes, therefore Ao = for all pairs. As t increases Ao becomes negative 
for some of the pairs of the bath modes, meaning that these modes become entangled. As 
time progresses, the number of entangled pairs saturates for all three couplings (Cf Fig. |Hl 
lower panel). Therefore, the increase in the absolute value of Ao is related primarily to the 
growth in population of the entangled modes. The maximum in the number of entangled 
modes for an early time may be associated with the creation of higher-order entanglement 
terms where three or more simultaneous excitations become important. Such higher-order 
entanglement is not captured by the Peres-Horodecki parameter. 

To characterize the degree of entanglement we use an additional measure, the entangle- 
ment of formation introduced by Wootters et al. j^, 48, 4^| (Cf. Appendix |Bj). For any 2® 2 
mixed state this quantity varies from zero (separable states) to one (maximally entangled 
states). Our results (Cf. Fig. EJ) are obtained by averaging the entanglement of formation 
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E(pij) (for two bath modes i and j) over all possible pairs of modes. The dynamics of (E(p)) 
is similar to those of the partial transpose parameter. It should be noted that the growth 
of entanglement shown in Figs. IHE1 is exclusive to the spin bath. 



V. SUMMARY AND CONCLUSIONS 



The similarities and differences of the relaxation dynamics of a primary system coupled to 
a spin or to a harmonic bath have been analyzed. The study was facilitated by the ability to 
obtain converged numerical results for finite period in time. In both cases this task becomes 
possible by employing a large but finite number of bath modes and controlling the degree 
of correlation. 



A. The similarities 



For all cases studied the extremely short time dynamics was identical. This period repre- 
sents the inertial response of the bath and is characterized by a zero derivative of the energy 
at the initial time t = 0. This non-Markovian dynamical evolution, "the slippage", is quite 
short and in many cases it can be ignored. The initial dynamics is closely related to the 
issue of the choice of the initial state. Preferably it should represent equilibrium system 
bath correlation and not be a product state. The Surrogate Hamiltonian method allows to 
create such a fully correlated initial state of the system-bath entity. However, in the present 
model differences in the dynamics between correlated and uncorrelated states seem to be 
insignificant, even in the strong coupling case. It is expected that the same phenomena 
would be observed in the harmonic bath. 

For weak system-bath coupling the dynamics induced by both baths are also similar. 
This is a numerical confirmation that in the weak coupling limit the harmonic bath can be 
mapped to the spin bath In addition in this limit the spin bath converges with 

only single excitations of the bath modes meaning that the system and bath are almost 
disentangled. This fact is consistent with the convergence to the Markovian limit jsfll ]. 

The decoherence properties of the harmonic and spin baths as determined by the loss of 
phase of cat states are found to be quite similar. This result is somewhat surprising since 
the ergodic properties of the two baths are different. To rationalize, one should notice that 
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the coherence in cat states composed of a superposition of two coherent states in a single 
mode does not represent entanglement. Therefore, this phase loss does not characterize 
decoherence in accordance with Alicki's notion jjyj. Moreover when a pure dephasing term 
was added it was found that it did not erode the phase coherence between cat states j^jj. 
We conclude that the decoherence properties of the two baths still require a further study. 



B. The differences 



The spin and harmonic baths begin to deviate when the initial excitation of the primary 
system is increased. This difference is observed for excitations where the dynamics generated 
by the harmonic bath is still Markovian. The first indication of differences is the necessity 
of two simultaneous bath excitation to converge the spin bath. For larger time periods, the 
spin bath saturates limiting the ability to assimilate the system's energy. The conclusion 
is that the limit of weak coupling is more restrictive in the spin bath case. The effect of 
saturation can be reduced if the bath Hamiltonian includes a mode-mode coupling term. 
This term causes diffusion of excitation between the modes, spreading the excitation over a 
greater number of bath modes. Thus bath modes, which are relatively far from resonance 
with the primary system become populated and the saturation is suppressed. Practically, 
this allows to increase the convergence timescale of the spin bath. 

In the medium coupling regime there is an overall good agreement between the two 
models. The spin bath, however, causes stronger relaxation, a fact, which becomes even 
more visible in the case of strong coupling. In this regime the deviations between the two 
baths become significant. 

The possibility of entanglement of bath modes mediated by the primary system is a major 
conceptual difference between the spin and harmonic baths. In the spin bath after a short 
initial period where only single excitations are excited, entanglement between pairs of spin 
sets in with what seems as an exponential growth. At later times the pair entanglement is 
replaced by higher order terms and the pair entanglement saturates. All these correlations 
are absent from the harmonic bath nevertheless the dynamics of the the primary systems are 
not very different except for extremely strong coupling. It could be possible that the in the 
present model the primary system and the system bath coupling terms are oversimplified. In 
addition the simulations should be extended to finite temperatures, where different behavior 
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of harmonic and spin baths is expected except for the weak coupling regime [5|, |3, |52j • The 
recently proposed random phase method |53] may be used to extend the above models to 
finite temperature applications. 

The present study is an important step in establishing the Surrogate Hamiltonian method 
as a practical simulation tool. The elucidation of the system-bath dynamics allows to tailor 
a simulation package in particular for ultrafast dynamical processes. 
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APPENDIX A: HARMONIC BATH VS SPIN BATH. 

The differences between harmonic and spin baths can be illuminated by studying the 
simple system of a spin-| coupled to a bath. For the harmonic bath the total Hamiltonian 
in second quantization is given by: 

H = hn&, + + V 6 ; + £})(*+ + *-) • (Al) 

3 3 

Thus the Heisenberg equations of motion for the system operators are: 

ja ± = 1[* ± , H] = ±itla ± - ia z J2 Ai(bi + b]) . (A2) 

3 

Similarly, the equations of motion for the bath operators are: 

— hj = —iujjhj — i\j(a + + cr_) . (A3) 

(Jib 

Since the annihilation and creation operators hj and b], satisfy the standard Bose commuta- 
tion relation, \bj, hj,] = 5jji, a closed set of equations is obtained for each of the independent 
bath modes. 

Now, let us consider a different model, where the primary system is coupled to a bath of 
spins- 1 . The total Hamiltonian may be written in the form: 

H = Oaj + ~ £ Uj&i - ~ + h.c. , (A4) 

3 3 
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where a x ,y,z designates the set of Pauli operators and d- l ± = cr % x ± cr % y are the usual ladder 
operators. For simplicity, we consider a system (0) consisting of a spin-| which interacts 
with a pair (1,2) of spins. Thus, the Hamiltonian of the whole system reads: 

H = 08° + i(u**J + co 2 al) - ^(d**i + di*l) - ^(alcrl + a <x 2 + ) . (A5) 

The Heisenberg equations of motion for the system operator are: 

= ±i(2Qai + ±d°,*l + ±d°dl) , (A6) 

and the equation of motion for the bath operator reads: 

d A.1,2 ,./ -1,2 . ^1,2 o s.l,2\ / A7 \ 

j t °± = ±i(a>i, 2 0± + — (T ± (T z ) ■ (A7) 

The commutation relations for spin operators are different from those of bosons: 

[a-i,&j] = -2ia k , (A8) 

which makes the set of the equations above non-closed. After some algebra, the equation 
for the bath operators becomes: 

« /-0'i.l,2\ ,. -0 - 1,2 , .-^1 -0 i -\ -1,2-2,1 -0-1, 2-2, In / a n \ 

-J^z^i ) = ±tMip(T g O-j. ±1 — 0± +1\2{V + (T± V- ~ <T-<T± ) • (A9) 

The triple correlations of the type a Q !r d- 1 ± a- 2 _ is a manifestation of the build-up of quantum 
entanglement - a specific correlation between different modes, which has no analogy in clas- 
sical physics. These correlations make a difference between the spin bath and the harmonic 
oscillator one, since the latter does not have quantum correlations between the bath modes. 



APPENDIX B: ENTANGLEMENT BETWEEN DIFFERENT BATH MODES. 

In order to check whether the reduced two-system density matrix p is entangled, first 
we will use the partial transposition criterion proposed by Peres j3] and Horodecki et al. 

A mixed state described by density matrix p is non-separable (and therefore cannot 
be written as a product state of two subsystems (i) and (j), p = p { (g) p-), iff the partial 
transposition of p with respect to one of the two subsystems has negative eigenvalues. The 
partial transpose p Tj is obtained by transposing in a matrix representation of p only those 
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indices corresponding to subsystem (j), i.e. pm^nu = Pmvw The following notation for 
matrix elements of a composite system is used: 

Pm M ,n. = (e m ® / M |p|e n <g> /„) , (Bl) 

where e m and / M denote the arbitrary orthonormal bases in Hilbert space describing the first 
(i) and second (j) system, respectively. 

Checking the positivity of the partial transpose is equivalent to checking the signs of the 
eigenvalues of p Tj or , alternatively, the signs of the following determinants: 

Wl = Pll, 11^22,22 — ^11,22^22,11) ^2 = Pl2, VlP%\, 21 — Pl2,2lP21,12 • (B2) 

In the case when one of the above determinants is negative, the state p is non-separable, 
and hence there is entanglement between the two subsystems p h p y 

Another entanglement measure for a mixed state of two spin-| particles is the entangle- 
ment of formation 42, 3,13- Explicitly, for the reduced two-system density matrix p the 



entanglement of formation is defined by 



E{p) = h I - [1 + y/1 - C(p) 2 \ 1 , (B3) 

where h is the binary entropy function h(x) = —xlog 2 x — (1 — x) log 2 (l — x) and C(p) is 
the concurrence. The concurrence is calculated in the following way: first we define the 
"spin-flipped" density matrix to be 

p=(a y ® a y )p*{a y ® Oy) , (B4) 

where the asterisk denotes complex conjugation of p in the standard basis 
{|00), 1 01) , 1 10) , 1 1 1) , } and a y expressed in the same basis is the matrix 

(o -A 

"y = • B5 

V / 

As both p and p are positive operators, it follows that the product pp, though non-Hermitian, 
also has only real and non-negative eigenvalues. Let the square roots of these eigenvalues, 
in decreasing order , be Ax, A 2 , A 3 , and A 4 . Then the concurrence of the density matrix p is 
defined as C = max{Ai — A 2 — A3 — A4, 0}. It should be noted that C = corresponds to an 
unentangled state, while C = 1 - to a completely entangled state and the entanglement of 
formation E is a monotonically increasing function of C. 
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